flowchart A[Define DAG] --> B(24 hours previous rain) A --> C(48 hours previous rain) B --> D(Define adjustment sets for each predictor) C --> D D --> E(Run all models satisfying\nthe back door criterion) E --> F(Average posterior probabilities) F --> G(Combine models in a single graph) style A fill:#44015466 style B fill:#3E4A894D style C fill:#3E4A894D style D fill:#26828E4D style E fill:#6DCD594D style F fill:#FDE7254D style G fill:#31688E4D
Statistical analysis using causal inference
Agalychnis lemur
Purpose
- Determine the adjustment sets that allow to infer a causal effect of environmental variables on vocal activity
- Evaluate causal effect of environmental factors on vocal activity of A. lemur with bayesian regression models
Analysis flowchart
Load packages
Code
# Unload packages
out <- sapply(paste('package:', names(sessionInfo()$otherPkgs), sep = ""), function(x) try(detach(x, unload = FALSE, character.only = TRUE), silent = T))
pkgs <-
c(
"remotes",
"viridis",
"brms",
"cowplot",
"posterior",
"readxl",
"HDInterval",
"kableExtra",
"knitr",
"ggdist",
"ggplot2",
"lunar",
"cowplot",
"maRce10/brmsish",
"warbleR",
"ohun",
"dagitty",
"ggdag",
"tidybayes",
"pbapply"
)
# install/ load packages
sketchy::load_packages(packages = pkgs)
options("digits" = 6, "digits.secs" = 5, knitr.table.format = "html")
# set evaluation false
opts_chunk$set(fig.width = 10, fig.height = 6, warning = FALSE, message = FALSE, tidy = TRUE)
# set working directory as project directory or one directory above,
opts_knit$set(root.dir = "..")
source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
source("~/Dropbox/R_package_testing/brmsish/R/check_rds_fits.R")
source("~/Dropbox/R_package_testing/brmsish/R/draw_extended_summary.R")Custom functions
Code
print_data_frame <- function(x){
# print table as kable
kb <-kable(x, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)
}
adjustment_set_formulas <-
function(dag,
exposure,
required_variable,
outcome,
effect = "total",
type = "minimal",
formula_parts = c(outcome),
latent = NULL,
remove = NULL,
plot = TRUE,
...) {
if (plot)
gg <- ggdag_adjustment_set(
.tdy_dag = tidy_dagitty(dag),
exposure = exposure,
outcome = outcome,
...
) + theme_dag()
temp_set <-
adjustmentSets(
x = dag,
exposure = exposure,
outcome = outcome,
effect = effect,
type = type
)
form_set <- lapply(temp_set, function(x) {
if (!is.null(remove))
x <- x[!x %in% remove]
form <-
paste(
formula_parts[1],
" ~ ",
exposure,
" + ",
paste(x, collapse = " + "),
if (length(formula_parts) == 2)
formula_parts[2]
else
NULL
)
return(form)
})
form_set <- form_set[!duplicated(form_set)]
if (!is.null(latent))
for (i in latent)
form_set <- form_set[!grepl(paste0(" ", i, " "), form_set)]
# form_set <- sapply(form_set, as.formula)
names(form_set) <- seq_along(form_set)
# add formula as attribute
attributes(form_set)$exposure.formula <- paste(formula_parts[1],
" ~ ",
exposure,
if (length(formula_parts) == 2)
formula_parts[2]
else
NULL)
if (plot)
return(gg) else
return(form_set)
}
# Define a function to remove special characters
remove_special_chars <- function(text) {
# Replace special characters with hyphen
cleaned_text <- gsub("[^a-zA-Z0-9]+", "-", text)
# Remove leading and trailing hyphens
cleaned_text <- gsub("^-+|-+$", "", cleaned_text)
return(cleaned_text)
}
pa <- function(...)
brms::posterior_average(...)
# to get average posterior values from models with different formulas
averaged_model <-
function(formulas,
data,
model_call,
ndraws = 1000,
save.files = TRUE,
path = ".",
suffix = NULL,
cores = 1,
name = NULL) {
# if (dir.exists(file.path(path, name))){
# cat("Directory already existed. Attempting to fit missing models\n")
# cat("Fitting models (step 1 out of 2) ...")
# } else
# dir.create(path = file.path(path, name))
cat("Fitting models (step 1 out of 2) ...")
fit_list <-
pblapply_brmsish_int(X = formulas, cl = cores, function(y) {
# make file name without special characters
mod_name <-
paste0(remove_special_chars(as.character(y)), ".RDS")
if (save.files &
!file.exists(file.path(path, mod_name))) {
cat("Fitting", y, "\n")
mc <-
gsub(pattern = "formula",
replacement = as.character(y),
x = model_call)
mc <- parse(text = mc)
fit <- eval(mc)
if (save.files)
saveRDS(fit, file = file.path(path, mod_name))
} else {
cat("Reading", y, "(already existed)\n")
fit <- readRDS(file.path(path, mod_name))
}
return(fit)
})
if (length(formulas) > 1){
cat("Averaging models (step 2 out of 2) ...")
average_call <-
parse(text = paste("pa(",
paste(
paste0("fit_list[[", seq_along(fit_list), "]]"), collapse = ", "
),
", ndraws = ",
ndraws,
")"))
# Evaluate the expression to create the call object
average_eval <- eval(average_call)
# add formula as attribute
attributes(average_eval)$averaged_fit_formulas <- formulas
rds_name <- if (is.null(suffix)) file.path(path, paste0(name, ".RDS")) else
file.path(path, paste0(suffix, "_", name, ".RDS"))
if (save.files)
saveRDS(average_eval, file = rds_name)
# return draws from average models
return(average_eval)
} else
cat("No model averaging conducted as a single formula was supplied")
}
to_change_percentage <- function(x){
(exp(x) - 1) * 100
}1 Prepare data
1.1 Read data
Code
clim_dat_2020 <- read_excel("./data/datos_metereologicos/Estacion_met/Dat_met_estacion_2022-01-11.xlsx",sheet = "2020")
clim_dat_2019 <- read_excel("./data/datos_metereologicos/Estacion_met/Dat_met_estacion_2022-01-11.xlsx",sheet = "2019")
clim_dat <- rbind(clim_dat_2019, clim_dat_2020)
clim_dat <- clim_dat[, c("filename", "Año", "Mes", "Día", "Hora", "Temp (°C)", "Humedad Relat.", "Precipitación")]
names(clim_dat) <- c("filename", "year", "month", "day", "hour", "temp", "HR", "rain")
clim_dat <- aggregate(cbind(rain, temp, HR) ~ filename + year + month + day + hour, clim_dat, mean)
clim_dat$year <- clim_dat$year + 2000
clim_dat$date <- as.Date(paste(clim_dat$year, clim_dat$month, clim_dat$day, sep = "-"))
clim_dat$date_hour <- paste(gsub("-", "", clim_dat$date), clim_dat$hour, sep = "-")
call_dat <- read.csv("./data/processed/call_rate_per_date_time_and_site.csv")
# remove 5 pm data and keep only form Sukia
call_dat <- call_dat[call_dat$hour != 17, ]
call_dat <- call_dat[call_dat$site != "LAGCHIMU", ]
call_dat$date_hour <- paste(sapply(as.character(call_dat$site_date_hour), function(x) strsplit(x, "_")[[1]][2]), call_dat$hour, sep = "-")
call_dat$temp <- sapply(1:nrow(call_dat), function(x){
y <- clim_dat$temp[clim_dat$date_hour == call_dat$date_hour[x]]
if (length(y) < 1) y <- NA
return(y)
}
)
call_dat$HR <- sapply(1:nrow(call_dat), function(x){
y <- clim_dat$HR[clim_dat$date_hour == call_dat$date_hour[x]]
if (length(y) < 1) y <- NA
return(y)
}
)
call_dat$rain <- sapply(1:nrow(call_dat), function(x){
y <- clim_dat$rain[clim_dat$date_hour == call_dat$date_hour[x]]
if (length(y) < 1) y <- NA
return(y)
}
)
# proportion of acoustic data with climatic data
sum(call_dat$date_hour %in% clim_dat$date_hour) / nrow(call_dat)
sum(!is.na(call_dat$temp)) / nrow(call_dat)
call_dat <- call_dat[!is.na(call_dat$temp), ]
call_dat$day <- as.numeric(substr(sapply(as.character(call_dat$site_date_hour), function(x) strsplit(x, "_")[[1]][2]), 7, 8))
call_dat$date <- as.Date(paste(call_dat$year, call_dat$month, call_dat$day, sep = "-"))
call_dat$moon.date <- ifelse(call_dat$hour < 12, as.Date(call_dat$date - 1), as.Date(call_dat$date))
call_dat$moon.date <- as.Date(call_dat$moon.date, origin = "1970-01-02")
## add moon
call_dat$moonlight <- lunar.illumination(call_dat$moon.date, shift = -6)
call_dat$date_hour_min <- strptime(paste(paste(call_dat$year, call_dat$month, call_dat$day, sep = "-"), paste(call_dat$hour, "00", sep = ":")), format="%Y-%m-%d %H:%M")
call_dat$hour_diff <- as.numeric(call_dat$date_hour_min - min(call_dat$date_hour_min)) / 3600
call_dat$rain_24 <- sapply(1:nrow(call_dat), function(x) sum(clim_dat$rain[strptime(clim_dat$date, format="%Y-%m-%d") == (strptime(call_dat$date[x], format="%Y-%m-%d") - 60 * 60 * 24)]))
call_dat$rain_48 <- sapply(1:nrow(call_dat), function(x) sum(clim_dat$rain[strptime(clim_dat$date, format="%Y-%m-%d") == (strptime(call_dat$date[x], format="%Y-%m-%d") - 60 * 60 * 48)]))
clim_dat$date_hour_min <- strptime(paste(paste(clim_dat$year, clim_dat$month, clim_dat$day, sep = "-"), paste(clim_dat$hour, "00", sep = ":")), format="%Y-%m-%d %H:%M")
clim_dat$hour_diff <- as.numeric(clim_dat$date_hour_min - min(call_dat$date_hour_min)) / 3600
# call_dat$prev_temp <- sapply(1:nrow(call_dat), function(x){
#
# # if(call_dat$hour_diff[x] < 48)
# # pt <- NA else
# pt <- mean(clim_dat$temp[clim_dat$hour_diff %in% (call_dat$hour_diff[x] - 48):(call_dat$hour_diff[x] - 24)])
#
# return(pt)
# })
write.csv(call_dat, "./data/processed/acoustic_and_climatic_data_by_hour.csv")2 Descriptive stats
Code
call_dat_site <- read.csv("./data/processed/call_rate_per_date_time_and_site.csv")- Total number of recordings: 9681
- Total recordings per site:
Code
agg_recs <- aggregate(rec_time ~ site, data = call_dat_site, length)
names(agg_recs)[1:2] <- c("site", "rec_count")
# print table as kable
kb <-kable(agg_recs, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)| site | rec_count | |
|---|---|---|
| 1 | LAGCHIMU | 5127 |
| 2 | SUKIA | 4554 |
Total recording time: 3224
Total recording time per site:
Code
agg_recs <- aggregate(rec_time ~ site, data = call_dat_site, sum)
names(agg_recs)[1:2] <- c("site", "recording_time")
agg_recs$recording_time <- round((agg_recs$recording_time) / 60)
# print table as kable
kb <-kable(agg_recs, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)| site | recording_time | |
|---|---|---|
| 1 | LAGCHIMU | 1707 |
| 2 | SUKIA | 1517 |
Total detections: 540203
Total detections per site:
Code
agg_recs <- aggregate(n_call ~ site, data = call_dat_site, sum)
names(agg_recs)[1:2] <- c("calls", "count")
# print table as kable
kb <-kable(agg_recs, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)| calls | count | |
|---|---|---|
| 1 | LAGCHIMU | 169485 |
| 2 | SUKIA | 370718 |
Call rate: 18.598587
Call rate per site:
Code
agg_recs <- aggregate(call_rate ~ site, data = call_dat_site, mean)
agg_recs$sd <- aggregate(call_rate ~ site, data = call_dat_site, FUN = stats::sd)[,2]
names(agg_recs)[1:3] <- c("site", "call_rate", "sd")
print_data_frame(agg_recs)| site | call_rate | sd | |
|---|---|---|---|
| 1 | LAGCHIMU | 11.018 | 18.039 |
| 2 | SUKIA | 27.133 | 87.778 |
Code
call_rate_hour <- read.csv("./data/processed/acoustic_and_climatic_data_by_hour.csv")
agg <- aggregate(cbind(temp, prev_temp, HR, rain, rain_24, rain_48, moonlight)~1, call_rate_hour, function(x) round(c(mean(x), sd(x), min(x), max(x)), 3))
agg <- as.data.frame(matrix(unlist(agg), ncol = 4, byrow = TRUE, dimnames = list(c("Temperature", "Previous temperature", "Relative humidity", "Night rain", "Rain 24 hours", "Rain 48 hours", "Moonlight"),c("mean", "sd", "min", "max"))))
# print table as kable
kb <-kable(agg, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)| mean | sd | min | max | |
|---|---|---|---|---|
| Temperature | 24.229 | 2.126 | 19.472 | 31.763 |
| Previous temperature | 24.097 | 0.988 | 20.005 | 27.111 |
| Relative humidity | 90.056 | 8.464 | 55.158 | 99.940 |
| Night rain | 0.079 | 0.457 | 0.000 | 10.417 |
| Rain 24 hours | 1.441 | 3.161 | 0.000 | 25.620 |
| Rain 48 hours | 1.453 | 3.170 | 0.000 | 25.620 |
| Moonlight | 0.495 | 0.355 | 0.000 | 1.000 |
Mean and SD temperature: 24.229179 (2.126374)
Mean previous temperature: 24.097073
Mean temperature: 24.229179
Mean cumulative rain per hour: 0.079232
Mean cumulative rain per hour previous 24 hours: 1.44064
Mean daily cumulative rain per hour previous 48 hours: 1.45317
3 Directed acyclical graphs (DAGs)
Code
coords <- list(
x = c(
sc_rain = -0.4,
evotranspiration = 0.5,
sc_prev_rain = 0.7,
sc_temp = -0.8,
sc_HR = 0,
n_call = 0,
sc_moonlight = 0.3,
hour = -0.5
),
y = c(
sc_rain = 0.4,
evotranspiration = 0.3,
sc_prev_rain = -0.5,
sc_temp = 0,
climate = 1,
sc_HR = -0.6,
n_call = 0,
sc_moonlight = 1,
hour = 0.9
)
)
# sc_temp + sc_HR + sc_moonlight + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour
# sc_temp = temp y meanT = prev_temp
dag_l <- dagify(sc_rain ~ evotranspiration,
sc_prev_rain ~ evotranspiration,
sc_temp ~ climate,
sc_temp ~ sc_rain,
sc_HR ~ sc_rain,
n_call ~ sc_HR,
n_call ~ hour,
n_call ~ sc_moonlight,
sc_moonlight ~ hour,
sc_temp ~ hour,
sc_HR ~ sc_temp,
sc_HR ~ sc_prev_rain,
sc_HR ~ sc_rain,
n_call ~ sc_temp,
n_call ~ sc_prev_rain,
n_call ~ sc_rain,
labels = c("n_call" = "Call\nrate", "sc_HR" = "Relative\nhumidity","sc_rain" = "Current\nrain", "sc_prev_rain" = "Prior\nrain", "sc_moonlight" = "Moonlight", "hour" = "Earth\nrotation", "sc_temp" = "Tempera-\nture", "evotranspiration" = "Evotrans-\npiration", "climate" = "Climate", latent = c("evotranspiration", "climate"), outcome = "n_call"), coords = coords)
tidy_dag <- tidy_dagitty(dag_l)
tidy_dag$data$type <- ifelse(is.na(tidy_dag$data$to), "outcome", "predictor")
tidy_dag$data$type[tidy_dag$data$name %in% c("evotranspiration", "climate")] <- "latent"
dat <- tidy_dag$data
shorten_distance <- c(0.07, 0.07)
dat$slope <- (dat$yend - dat$y) / (dat$xend - dat$x)
distance <- sqrt((dat$xend - dat$x)^2 + (dat$yend - dat$y)^2)
proportion <- shorten_distance[1]/distance
dat$xend <- (1 - proportion/2) * dat$xend + (proportion/2 * dat$x)
dat$yend <- (1 - proportion/2) * dat$yend + (proportion/2*dat$y)
proportion <- shorten_distance[2]/distance
dat$xstart <- (1-proportion/2)*(dat$x - dat$xend) + dat$xend
dat$ystart <- (1-proportion/2)*(dat$y - dat$yend) + dat$yend
tidy_dag$data <- dat
basic_dag <- ggplot(tidy_dag, aes(
x = x,
y = y,
xend = xend,
yend = yend
)) +
scale_color_viridis_d(begin = 0.2, end = 0.8, alpha = 0.5) +
geom_dag_text(color = "black", aes(label = label, color = label)) + labs(color = "Type") +
theme_dag() + theme(legend.position = "bottom") + guides(colour = guide_legend(override.aes = list(size =
10))) + geom_dag_point(aes(color = type), size = 30) + expand_limits(y = c(-0.67, 1.1)) +
geom_dag_edges_fan(
edge_color = viridis(10, alpha = 0.4)[2],
arrow = grid::arrow(length = grid::unit(10, "pt"), type = "closed"),
aes(
x = xstart,
y = ystart,
xend = xend,
yend = yend
)
)
ggsave(plot = basic_dag, filename = "./output/dag_300dpi_v02.png", width = 9, height = 5, dpi = 300, bg = "white")
dag_24 <- dagify(sc_rain ~ evotranspiration,
sc_rain_24 ~ evotranspiration,
sc_temp ~ climate,
sc_temp ~ sc_rain,
sc_HR ~ sc_rain,
n_call ~ sc_HR,
n_call ~ hour,
n_call ~ sc_moonlight,
sc_moonlight ~ hour,
sc_temp ~ hour,
sc_HR ~ sc_temp,
sc_HR ~ sc_rain_24,
sc_HR ~ sc_rain,
n_call ~ sc_temp,
n_call ~ sc_rain_24,
n_call ~ sc_rain,
labels = c("n_call" = "Call rate", "sc_HR" = "Relative humidity","sc_rain" = "Night Rain", "sc_rain_24" = "Previous Rain", "sc_moonlight" = "Moonlight", "hour" = "Earth rotation", "sc_temp" = "Temperature", "evotranspiration" = "Evotranspiration", "climate" = "Climate", latent = c("evotranspiration", "climate"), outcome = "n_call"))
dag_48 <- dagify(sc_rain ~ evotranspiration,
sc_rain_48 ~ evotranspiration,
sc_temp ~ climate,
sc_temp ~ sc_rain,
sc_HR ~ sc_rain,
n_call ~ sc_HR,
n_call ~ hour,
n_call ~ sc_moonlight,
sc_moonlight ~ hour,
sc_temp ~ hour,
sc_HR ~ sc_temp,
sc_HR ~ sc_rain_48,
sc_HR ~ sc_rain,
n_call ~ sc_temp,
n_call ~ sc_rain_48,
n_call ~ sc_rain,
labels = c("n_call" = "Call rate", "sc_HR" = "Relative humidity","sc_rain" = "Night Rain", "sc_rain_48" = "Previous Rain", "sc_moonlight" = "Moonlight", "hour" = "Earth rotation", "sc_temp" = "Temperature", "evotranspiration" = "Evotranspiration", "climate" = "Climate", latent = c("evotranspiration", "climate"), outcome = "n_call"))4 Bayesian regression models
4.1 Scale variables and set model parameters
Code
call_rate_hour <- read.csv("./data/processed/acoustic_and_climatic_data_by_hour.csv")
# make hour a factor
call_rate_hour$hour <- factor(call_rate_hour$hour)
# scale and mean-center
call_rate_hour$sc_temp <- scale(call_rate_hour$temp)
call_rate_hour$sc_HR <- scale(call_rate_hour$HR)
call_rate_hour$sc_rain <- scale(call_rate_hour$rain)
call_rate_hour$sc_rain_24 <- scale(call_rate_hour$rain_24)
call_rate_hour$sc_rain_48 <- scale(call_rate_hour$rain_48)
call_rate_hour$sc_moonlight <- scale(call_rate_hour$moonlight)
priors <- c(prior(normal(0, 4), class = "b"))
chains <- 4
iter <- 100004.2 Fit all models
Code
param_grid <- expand.grid(
dag = c("dag_24", "dag_48"),
exposure = c("sc_temp", "sc_HR", "sc_moonlight", "sc_rain", "sc_rain_24", "sc_rain_48"),
effect = c("total", "direct"), stringsAsFactors = FALSE
)
param_grid$name <- apply(param_grid, 1, paste, collapse = "-")
# remove wrong dags for previous rain
param_grid <- param_grid[!(param_grid$dag == "dag_24" & param_grid$exposure == "sc_rain_48") & !(param_grid$dag == "dag_48" & param_grid$exposure == "sc_rain_24"), ]
adjustment_sets_list <- pblapply(seq_len(nrow(param_grid)), cl = 1, function(x){
forms <- adjustment_set_formulas(
dag = if (param_grid$dag[x] == "dag_24") dag_24 else dag_48,
type = if (param_grid$effect[x] == "total") "all" else "minimal",
exposure = param_grid$exposure[x],
outcome = "n_call",
effect = param_grid$effect[x],
required_variable = "hour",
formula_parts = c(
"n_call | resp_rate(rec_time)",
"+ ar(p = 2, time = hour_diff, gr = hour)"
),
latent = c("evotranspiration", "climate"),
remove = "hour",
plot = FALSE
)
return(forms)
})
names(adjustment_sets_list) <- param_grid$name
param_grid$model <-sapply(seq_len(nrow(param_grid)), function(x) {
if (param_grid$effect[x] == "direct")
adjustment_sets_list[[which(names(adjustment_sets_list) == param_grid$name[x])]] else
NA
})
param_grid$model <- unlist(
param_grid$model)
param_grid <- as.data.frame(param_grid)
param_grid$model[!is.na(param_grid$model)] <- remove_special_chars(param_grid$model[!is.na(param_grid$model)] )
param_grid$model <- c("total_effect_temperature_with_rain_24", "total_effect_temperature_with_rain_48", "total_effect_humidity_with_rain_24", "total_effect_humidity_with_rain_48", "total_effect_moon_with_rain_24", "total_effect_moon_with_rain_48", "total_effect_rain_with_rain_24", "total_effect_rain_with_rain_48", "total_effect_previous_rain_24", "total_effect_previous_rain_48", param_grid$model[!is.na(param_grid$model)])
param_grid$exposure.name <- param_grid$exposure
param_grid$exposure.name[grep("temp", param_grid$exposure.name)] <- "Temperature"
param_grid$exposure.name[grep("HR", param_grid$exposure.name)] <- "Relative humidity"
param_grid$exposure.name[grep("moon", param_grid$exposure.name)] <- "Moonlight"
param_grid$exposure.name[grep("rain$", param_grid$exposure.name)] <- "Current rain"
param_grid$exposure.name[grep("rain_24", param_grid$exposure.name)] <- "Previous rain (24h)"
param_grid$exposure.name[grep("rain_48", param_grid$exposure.name)] <- "Previous rain (48h)"
table(param_grid$exposure.name)
write.csv(x = param_grid, file = "./data/processed/direct_and_total_effect_model_data_frame.csv", row.names = FALSE)
direct_adjustment_sets_list <- adjustment_sets_list[grep("direct", names(adjustment_sets_list))]
for(i in seq_along(direct_adjustment_sets_list))
pa_comb_mod <-
averaged_model(
formulas = direct_adjustment_sets_list[[i]],
data = call_rate_hour,
suffix = "direct",
model_call = "brm(formula, data = call_rate_hour, iter = iter, chains = chains, cores = chains, family = negbinomial(), prior = priors, backend = 'cmdstanr')",
save.files = TRUE,
path = "./data/processed/averaged_models",
# name = "temperature_with_rain_24",
cores = 1
)
model_call = "brm(formula, data = call_rate_hour, iter = iter, chains = chains, cores = chains, family = negbinomial(), prior = priors, backend = 'cmdstanr')"
formulas <- unlist(direct_adjustment_sets_list)
path <- "./data/processed/single_models"
fit_list <- pblapply_brmsish_int(X = formulas, cl = 1, function(y) {
# make file name without special characters
mod_name <-
paste0(remove_special_chars(as.character(y)), ".RDS")
if (!file.exists(file.path(path, mod_name))) {
cat("Fitting", y, "\n")
mc <-
gsub(pattern = "formula",
replacement = as.character(y),
x = model_call)
mc <- parse(text = mc)
fit <- eval(mc)
if (save.files)
saveRDS(fit, file = file.path(path, mod_name))
}
})5 Results
Code
param_grid <- read.csv(file = "./data/processed/direct_and_total_effect_model_data_frame.csv")
param_grid$files <- file.path("./data/processed/averaged_models", paste0(param_grid$model, ".RDS"))
for(i in unique(param_grid$exposure.name)){
Y <- param_grid[param_grid$exposure.name == i, ]
cat(paste("\n##", i), "\n")
cat("\n### Direct effects\n")
for(e in which(Y$effect == "direct")){
if (grepl("24", Y$model[e]))
cat("\n#### 24 hour previous rain model:\n") else
cat("\n#### 48 hour previous rain model:\n")
extended_summary(
read.file = Y$files[e],
highlight = TRUE,
remove.intercepts = TRUE,
print.name = FALSE
)
cat("\n")
}
cat("\n### Total effect\n")
for(w in which(Y$effect == "total")){
if (grepl("24", Y$files[w]))
cat("\n#### 24 hour previous rain model:\n") else
cat("\n#### 48 hour previous rain model:\n")
draws <- readRDS(Y$files[w])
draw_extended_summary(
draws,
highlight = TRUE,
remove.intercepts = TRUE
)
cat("\n")
cat("\n##### Summary of single models:\n")
# print summary
print(readRDS(gsub("\\.RDS", "_fit_summary.RDS", Y$files[w])))
}
cat("\n")
}5.1 Temperature
5.1.1 Direct effects
5.1.1.1 24 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 13808 | 14519.5 | 1687235322 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_temp | 0.577 | 0.498 | 0.656 | 1 | 13808.0 | 14685.6 |
| b_sc_HR | 0.314 | 0.242 | 0.387 | 1 | 18992.0 | 14519.5 |
| b_sc_rain | 0.037 | 0.002 | 0.074 | 1 | 23964.4 | 15266.3 |
| b_sc_rain_24 | 0.096 | 0.057 | 0.135 | 1 | 23319.6 | 15914.6 |
5.1.1.2 48 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 9460.49 | 12432.4 | 293904519 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_temp | 0.565 | 0.487 | 0.643 | 1 | 9460.49 | 12432.4 |
| b_sc_HR | 0.307 | 0.235 | 0.379 | 1 | 11994.69 | 13312.6 |
| b_sc_rain | 0.029 | -0.007 | 0.066 | 1 | 17566.34 | 16050.0 |
| b_sc_rain_48 | 0.006 | -0.032 | 0.044 | 1 | 17289.44 | 16222.0 |
5.1.2 Total effect
5.1.2.1 24 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_temp | 0.308 | 0.258 | 0.361 | 1.043 | 44.190 | 893.839 |
| b_sc_rain | 0.041 | 0.004 | 0.078 | 1.026 | 152.336 | 830.415 |
5.1.2.1.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2542.43 | 5383.76 | 629626678 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2238.34 | 5028.64 | 2096770826 |
| 3 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2320.90 | 4619.08 | 1430954384 |
| 4 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2292.54 | 4473.95 | 1881604582 |
5.1.2.2 48 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_temp | 0.305 | 0.258 | 0.353 | 1.005 | 961.292 | 983.283 |
| b_sc_rain | 0.036 | -0.002 | 0.076 | 1 | 1113.522 | 969.564 |
5.1.2.2.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2529.79 | 5016.37 | 1471048663 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_moonlight + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2588.68 | 5844.63 | 412112948 |
| 3 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2622.63 | 5641.91 | 41093237 |
| 4 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2334.61 | 4741.38 | 366114719 |
5.2 Relative humidity
5.2.1 Direct effects
5.2.1.1 24 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 11039.8 | 13170.1 | 1554924087 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_HR | 0.314 | 0.240 | 0.388 | 1 | 13006.8 | 13679.9 |
| b_sc_rain | 0.037 | 0.002 | 0.075 | 1 | 18235.1 | 15400.0 |
| b_sc_rain_24 | 0.096 | 0.057 | 0.135 | 1 | 16966.8 | 15476.0 |
| b_sc_temp | 0.578 | 0.499 | 0.657 | 1 | 11039.8 | 13170.1 |
5.2.1.2 48 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 8448.99 | 11500.1 | 1931567299 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_HR | 0.308 | 0.234 | 0.380 | 1 | 10649.49 | 12433.6 |
| b_sc_rain | 0.029 | -0.008 | 0.067 | 1 | 16664.06 | 15531.5 |
| b_sc_rain_48 | 0.005 | -0.032 | 0.043 | 1 | 15745.61 | 14835.4 |
| b_sc_temp | 0.566 | 0.487 | 0.645 | 1 | 8448.99 | 11500.1 |
5.2.2 Total effect
5.2.2.1 24 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_HR | 0.309 | 0.240 | 0.383 | 1.006 | 942.064 | 903.519 |
| b_sc_rain | 0.038 | 0.002 | 0.076 | 1.003 | 995.090 | 912.254 |
| b_sc_rain_24 | 0.094 | 0.057 | 0.133 | 1.003 | 993.441 | 772.255 |
| b_sc_temp | 0.575 | 0.496 | 0.655 | 1.004 | 1022.397 | 861.998 |
5.2.2.1.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_HR + sc_moonlight + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2556.01 | 5127.35 | 1982612331 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2182.06 | 4328.43 | 1554924087 |
5.2.2.2 48 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_HR | 0.303 | 0.230 | 0.371 | 1.006 | 801.527 | 931.946 |
| b_sc_rain | 0.030 | -0.008 | 0.066 | 1.002 | 812.238 | 1071.444 |
| b_sc_rain_48 | 0.003 | -0.034 | 0.044 | 1.005 | 889.913 | 984.434 |
| b_sc_temp | 0.564 | 0.480 | 0.644 | 1.003 | 821.550 | 773.876 |
5.2.2.2.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_HR + sc_moonlight + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2326.55 | 4696.55 | 846435536 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2166.28 | 4627.63 | 1931567299 |
5.3 Moonlight
5.3.1 Direct effects
5.3.1.1 48 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 30283.5 | 15929.8 | 1624275037 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_moonlight | -0.129 | -0.193 | -0.064 | 1 | 30283.5 | 15929.8 |
5.3.1.2 48 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 30283.5 | 15929.8 | 1624275037 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_moonlight | -0.129 | -0.193 | -0.064 | 1 | 30283.5 | 15929.8 |
5.3.2 Total effect
5.3.2.1 24 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_moonlight | -0.101 | -0.17 | -0.039 | 1.04 | 40.665 | 750.448 |
5.3.2.1.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2463.32 | 4856.49 | 1624275037 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2613.63 | 5435.51 | 1331671454 |
| 3 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2426.99 | 4881.26 | 1714115616 |
| 4 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2184.13 | 4557.11 | 946500998 |
| 5 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2386.65 | 5083.35 | 2088023938 |
| 6 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2571.23 | 4971.80 | 1587538237 |
| 7 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2319.91 | 4626.15 | 1236984196 |
| 8 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2197.61 | 4746.52 | 1987527846 |
| 9 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2353.89 | 5263.44 | 924618251 |
| 10 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2409.47 | 4652.12 | 1063522846 |
| 11 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2407.64 | 4851.64 | 259660598 |
| 12 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2473.17 | 4731.37 | 1625659188 |
| 13 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2596.10 | 5292.48 | 1259975215 |
| 14 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 1998.35 | 4047.25 | 192228769 |
| 15 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2365.76 | 4743.96 | 559459570 |
| 16 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2432.09 | 4364.34 | 1786523665 |
5.3.2.2 48 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_moonlight | -0.105 | -0.166 | -0.038 | 1.02 | 111.506 | 659.054 |
5.3.2.2.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + +ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2534.83 | 5170.40 | 169376128 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2664.91 | 5585.82 | 731739252 |
| 3 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2449.93 | 4958.79 | 230289957 |
| 4 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2283.15 | 4266.04 | 1735286727 |
| 5 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2486.75 | 4561.76 | 1539258099 |
| 6 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2525.61 | 5053.20 | 1234320719 |
| 7 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2481.99 | 4862.29 | 183236870 |
| 8 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2296.24 | 4114.83 | 239469808 |
| 9 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2298.87 | 4495.74 | 1126617762 |
| 10 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_HR + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2442.38 | 4594.40 | 684931243 |
| 11 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 1 (5e-05%) | 0 | 2416.00 | 4930.19 | 1093486403 |
| 12 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2146.30 | 4395.84 | 1670597608 |
| 13 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2450.19 | 4442.71 | 627408872 |
| 14 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2513.43 | 4849.68 | 1827601939 |
| 15 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2101.51 | 4415.14 | 104020182 |
| 16 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2495.11 | 5464.93 | 598459253 |
| 17 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + prev_temp + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2120.85 | 4487.97 | 2041361978 |
| 18 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2336.51 | 4816.06 | 1759350038 |
| 19 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 1 (5e-05%) | 0 | 2453.12 | 5040.07 | 927628753 |
| 20 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2360.45 | 5250.70 | 2028572807 |
| 21 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2642.43 | 5077.10 | 293339431 |
| 22 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2342.88 | 4887.79 | 1868346384 |
| 23 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2515.33 | 5176.64 | 552665161 |
| 24 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2253.53 | 4141.73 | 1764642823 |
| 25 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_HR + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2634.31 | 5119.83 | 505424120 |
| 26 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2427.59 | 5073.57 | 1761124445 |
| 27 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2609.57 | 5134.76 | 1073467634 |
| 28 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2461.43 | 4814.97 | 1419285187 |
| 29 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2596.36 | 5583.40 | 1232772932 |
| 30 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2271.50 | 5053.09 | 278133969 |
| 31 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 1 (5e-05%) | 0 | 2507.89 | 4315.91 | 866725413 |
| 32 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_moonlight + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2126.46 | 4549.18 | 1885323046 |
5.4 Current rain
5.4.1 Direct effects
5.4.1.1 24 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain + sc_HR + sc_rain_24 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 10228.7 | 12456.7 | 1281527469 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain | 0.037 | 0.001 | 0.075 | 1 | 18004.9 | 15909.5 |
| b_sc_HR | 0.313 | 0.240 | 0.386 | 1 | 13025.4 | 13472.6 |
| b_sc_rain_24 | 0.096 | 0.059 | 0.134 | 1 | 17444.9 | 15779.0 |
| b_sc_temp | 0.577 | 0.499 | 0.655 | 1 | 10228.7 | 12456.7 |
5.4.1.2 48 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain + sc_HR + sc_rain_48 + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 7371.36 | 11882.6 | 1129340659 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain | 0.029 | -0.007 | 0.066 | 1 | 14353.91 | 14479.2 |
| b_sc_HR | 0.307 | 0.235 | 0.380 | 1 | 8890.14 | 12279.5 |
| b_sc_rain_48 | 0.005 | -0.032 | 0.043 | 1 | 13843.68 | 14097.9 |
| b_sc_temp | 0.566 | 0.487 | 0.644 | 1 | 7371.36 | 11882.6 |
5.4.2 Total effect
5.4.2.1 24 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain | 0.000 | -0.034 | 0.035 | 1.001 | 942.79 | 942.018 |
| b_sc_rain_24 | 0.076 | 0.039 | 0.113 | 1 | 1056.52 | 982.799 |
5.4.2.1.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain + sc_moonlight + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2531.75 | 4835.39 | 1596946440 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain + sc_rain_24 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2369.16 | 4574.59 | 1687235322 |
5.4.2.2 48 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain | -0.005 | -0.042 | 0.030 | 1 | 1024.455 | 819.064 |
| b_sc_rain_48 | 0.011 | -0.027 | 0.049 | 1.009 | 932.089 | 975.233 |
5.4.2.2.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain + sc_moonlight + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2647.66 | 5129.91 | 1485832077 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain + sc_rain_48 + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2675.26 | 5099.07 | 167092432 |
5.5 Previous rain (24h)
5.5.1 Direct effects
5.5.1.1 24 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 11374.2 | 13445.2 | 1379661434 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain_24 | 0.096 | 0.058 | 0.135 | 1 | 20345.7 | 15663.3 |
| b_sc_HR | 0.314 | 0.240 | 0.386 | 1 | 14938.4 | 13445.2 |
| b_sc_rain | 0.038 | 0.001 | 0.074 | 1 | 21502.1 | 16005.3 |
| b_sc_temp | 0.578 | 0.498 | 0.656 | 1 | 11374.2 | 13493.9 |
5.5.2 Total effect
5.5.2.1 24 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain_24 | 0.088 | 0.051 | 0.127 | 1.01 | 883.631 | 733.445 |
| b_sc_rain | 0.043 | 0.005 | 0.081 | 1.001 | 908.203 | 446.195 |
5.5.2.1.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2692.30 | 5204.02 | 423216939 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2556.12 | 4635.90 | 610868624 |
| 3 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 1 (5e-05%) | 0 | 2559.09 | 4449.50 | 1988187279 |
| 4 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_24 + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2077.52 | 4134.53 | 640249946 |
5.6 Previous rain (48h)
5.6.1 Direct effects
5.6.1.1 48 hour previous rain model:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-inv_gamma(0.4, 0.3) | n_call | resp_rate(rec_time) ~ sc_rain_48 + sc_HR + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 11112 | 13056.2 | 1238572765 |
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain_48 | 0.006 | -0.031 | 0.043 | 1 | 17211.0 | 15210.0 |
| b_sc_HR | 0.308 | 0.236 | 0.380 | 1 | 14154.7 | 13956.7 |
| b_sc_rain | 0.029 | -0.007 | 0.066 | 1 | 19170.2 | 15964.2 |
| b_sc_temp | 0.566 | 0.486 | 0.645 | 1 | 11112.0 | 13056.2 |
5.6.2 Total effect
5.6.2.1 48 hour previous rain model:
| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| b_sc_rain_48 | 0.005 | -0.033 | 0.043 | 1.014 | 711.511 | 873.701 |
| b_sc_rain | 0.036 | 0.001 | 0.073 | 1 | 857.372 | 792.928 |
5.6.2.1.1 Summary of single models:
| priors | formula | iterations | chains | thinning | warmup | diverg_transitions | rhats > 1.05 | min_bulk_ESS | min_tail_ESS | seed | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_moonlight + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2447.89 | 5823.08 | 1081516016 |
| 2 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 2 (1e-04%) | 0 | 2413.93 | 4851.86 | 217409394 |
| 3 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_rain + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2506.36 | 5173.40 | 724206553 |
| 4 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_48 + prev_temp + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2641.41 | 5493.32 | 509251240 |
| 5 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_48 + sc_moonlight + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2386.41 | 4640.41 | 518390312 |
| 6 | b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0, 2.5) shape-gamma(0.01, 0.01) | n_call | resp_rate(rec_time) ~ sc_rain_48 + sc_rain + sc_temp + ar(p = 2, time = hour_diff, gr = hour) | 10000 | 4 | 1 | 5000 | 0 (0%) | 0 | 2575.04 | 5332.87 | 76427694 |
5.7 Combined results with causal inference estimates
Takes the posterior probability distributions from the right causal models
Code
param_grid <- read.csv(file = "./data/processed/direct_and_total_effect_model_data_frame.csv")
param_grid <- param_grid[param_grid$effect == "direct",]
param_grid$file <- paste0(remove_special_chars(param_grid$model), ".RDS")
rdss_24 <- list.files("./data/processed/averaged_models", pattern = "24.RDS", full.names = TRUE)
combined_draws_list <- lapply(rdss_24, function(x) {
total_draws <- readRDS(x)
exp <-
attributes((attributes(total_draws)$averaged_fit_formulas))$exposure.formula
exp <-
gsub("n_call | resp_rate(rec_time) ~ ", "", exp, fixed = TRUE)
exposure <-
exp <-
gsub(" + ar(p = 2, time = hour_diff, gr = hour)", "", exp, fixed = TRUE)
exp <- paste0("b_", exp)
total_draws <-
total_draws[, colnames(total_draws) == exp, drop = FALSE]
names(total_draws) <- exp
direct_fit_file <-
param_grid$file[param_grid$exposure == exposure]
direct_fit_file <- direct_fit_file[!duplicated(direct_fit_file)]
if (length(direct_fit_file) > 1)
direct_fit_file <- grep("24", direct_fit_file, value = TRUE)
direct_fit <-
readRDS(file = file.path("./data/processed/averaged_models", direct_fit_file))
direct_draws <-
posterior::merge_chains(as_draws(direct_fit, variable = exp))
direct_draws <-
as.data.frame(thin_draws(direct_draws, thin = length(direct_draws[[1]][[1]])
/ (nrow(total_draws)))[[1]])
direct_draws$effect <- "direct"
total_draws$effect <- "total"
draws <- rbind(direct_draws, total_draws)
return(draws)
})
combined_draws <- do.call(cbind, combined_draws_list)
combined_draws <- combined_draws[, c(which(sapply(combined_draws, is.numeric)), ncol(combined_draws))]
combined_draws[,-ncol(combined_draws)] <- to_change_percentage(combined_draws[,-ncol(combined_draws)])
# combined_draws <- as.data.frame(combined_draws)
combined_draws$effect <- ifelse(combined_draws$effect == "direct", "Direct", "Total")
saveRDS(combined_draws, "./data/processed/combined_draws_for_total_and_direct_effects_24h_previous_rain.RDS")Code
combined_draws <- readRDS("./data/processed/combined_draws_for_total_and_direct_effects_24h_previous_rain.RDS")
fill_colors <- viridis::mako(10)[c(8, 4)]
gg_dists <- draw_extended_summary(
draws = combined_draws,
highlight = TRUE,
remove.intercepts = TRUE,
fill = adjustcolor(fill_colors, alpha.f = 0.4),
by = "effect",
gsub.pattern = c(
"b_sc_HR",
"b_sc_rain$",
"b_sc_rain_24",
"b_sc_temp",
"b_sc_moonlight"
),
gsub.replacement = c(
"Relative\nhumidity",
"Current\nrain",
"Prior\nrain",
"Temperature",
"Moonlight"
),
ylab = "Variable",
xlab = "Effect size (% of change)"
)| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| Direct-Current rain | 3.652 | -0.002 | 7.942 | 1.009 | 999.221 | 949.102 |
| Total-Current rain | 0.000 | -3.373 | 3.596 | 1.001 | 942.790 | 942.018 |
| Direct-Moonlight | -12.018 | -17.074 | -6.217 | 1.001 | 972.135 | 944.341 |
| Total-Moonlight | -9.600 | -15.624 | -3.860 | 1.04 | 40.665 | 750.448 |
| Direct-Prior rain | 10.060 | 5.793 | 14.735 | 1 | 916.714 | 983.283 |
| Total-Prior rain | 9.189 | 5.189 | 13.581 | 1.01 | 883.631 | 733.445 |
| Direct-Relative humidity | 37.000 | 27.013 | 47.025 | 1.001 | 902.884 | 848.891 |
| Total-Relative humidity | 36.187 | 27.104 | 46.613 | 1.006 | 942.064 | 903.519 |
| Direct-Temperature | 77.591 | 65.145 | 92.509 | 1.002 | 955.952 | 875.689 |
| Total-Temperature | 36.112 | 29.497 | 43.428 | 1.043 | 44.190 | 893.839 |
Code
gg_dists + ggplot2::scale_fill_manual(values = fill_colors) +
ggplot2::theme(
axis.ticks.length = ggplot2::unit(0, "pt"),
plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"),
legend.position = "inside",
legend.position.inside = c(0.7, 0.7)
)Code
ggsave("./output/posterior_distribution_change_percentage_24h_previous_rain.png", width = 5, height = 3, dpi = 300)Code
param_grid <- read.csv(file = "./data/processed/direct_and_total_effect_model_data_frame.csv")
param_grid <- param_grid[param_grid$effect == "direct",]
param_grid$file <- paste0(remove_special_chars(param_grid$model), ".RDS")
rdss_48 <- list.files("./data/processed/averaged_models", pattern = "48.RDS", full.names = TRUE)
combined_draws_list <- lapply(rdss_48, function(x) {
total_draws <- readRDS(x)
exp <-
attributes((attributes(total_draws)$averaged_fit_formulas))$exposure.formula
exp <-
gsub("n_call | resp_rate(rec_time) ~ ", "", exp, fixed = TRUE)
exposure <-
exp <-
gsub(" + ar(p = 2, time = hour_diff, gr = hour)", "", exp, fixed = TRUE)
exp <- paste0("b_", exp)
total_draws <-
total_draws[, colnames(total_draws) == exp, drop = FALSE]
names(total_draws) <- exp
direct_fit_file <-
param_grid$file[param_grid$exposure == exposure]
direct_fit_file <- direct_fit_file[!duplicated(direct_fit_file)]
if (length(direct_fit_file) > 1)
direct_fit_file <- grep("48", direct_fit_file, value = TRUE)
direct_fit <-
readRDS(file = file.path("./data/processed/averaged_models", direct_fit_file))
direct_draws <-
posterior::merge_chains(as_draws(direct_fit, variable = exp))
direct_draws <-
as.data.frame(thin_draws(direct_draws, thin = length(direct_draws[[1]][[1]])
/ (nrow(total_draws)))[[1]])
direct_draws$effect <- "direct"
total_draws$effect <- "total"
draws <- rbind(direct_draws, total_draws)
return(draws)
})
combined_draws <- do.call(cbind, combined_draws_list)
combined_draws <- combined_draws[, c(which(sapply(combined_draws, is.numeric)), ncol(combined_draws))]
combined_draws[,-ncol(combined_draws)] <- to_change_percentage(combined_draws[,-ncol(combined_draws)])
# combined_draws <- as.data.frame(combined_draws)
combined_draws$effect <- ifelse(combined_draws$effect == "direct", "Direct", "Total")
saveRDS(combined_draws, "./data/processed/combined_draws_for_total_and_direct_effects_48h_previous_rain.RDS")Code
combined_draws <- readRDS("./data/processed/combined_draws_for_total_and_direct_effects_48h_previous_rain.RDS")
gg_dists <- draw_extended_summary(
draws = combined_draws,
highlight = TRUE,
remove.intercepts = TRUE,
fill = adjustcolor(fill_colors, alpha.f = 0.4),
by = "effect",
gsub.pattern = c(
"b_sc_HR",
"b_sc_rain$",
"b_sc_rain_48",
"b_sc_temp",
"b_sc_moonlight"
),
gsub.replacement = c(
"Relative\nhumidity",
"Current\nrain",
"Prior\nrain",
"Temperature",
"Moonlight"
),
ylab = "Variable",
xlab = "Effect size (% of change)"
)| Estimate | l-95% CI | u-95% CI | Rhat | Bulk_ESS | Tail_ESS | |
|---|---|---|---|---|---|---|
| Direct-Current rain | 2.940 | -0.477 | 6.832 | 1 | 933.292 | 847.036 |
| Total-Current rain | -0.544 | -4.093 | 3.071 | 1 | 1024.455 | 819.064 |
| Direct-Moonlight | -12.018 | -17.074 | -6.217 | 1.001 | 972.135 | 944.341 |
| Total-Moonlight | -9.925 | -15.276 | -3.765 | 1.02 | 111.506 | 659.054 |
| Direct-Prior rain | 0.623 | -2.979 | 4.124 | 1.003 | 857.680 | 1009.747 |
| Total-Prior rain | 0.540 | -3.268 | 4.427 | 1.014 | 711.511 | 873.701 |
| Direct-Relative humidity | 36.477 | 26.300 | 46.337 | 0.999 | 1043.482 | 878.337 |
| Total-Relative humidity | 35.405 | 25.805 | 44.967 | 1.006 | 801.527 | 931.946 |
| Direct-Temperature | 76.199 | 63.368 | 89.305 | 1.004 | 1009.050 | 835.503 |
| Total-Temperature | 35.603 | 29.387 | 42.312 | 1.005 | 961.292 | 983.283 |
Code
gg_dists + ggplot2::scale_fill_manual(values = fill_colors) +
ggplot2::theme(
axis.ticks.length = ggplot2::unit(0, "pt"),
plot.margin = ggplot2::margin(0, 0, 0, 0, "pt"),
legend.position = "inside",
legend.position.inside = c(0.7, 0.7)
)Code
ggsave("./output/posterior_distribution_change_percentage_48h_previous_rain.png", width = 5, height = 3, dpi = 300)Takeaways
- Variation in call activity strongly linked to environmental variation
Session information
R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=es_CR.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
time zone: America/Costa_Rica
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] pbapply_1.7-2 tidybayes_3.0.6 ggdag_0.2.12 dagitty_0.3-4
[5] ohun_1.0.1 warbleR_1.1.30 NatureSounds_1.0.4 seewave_2.2.3
[9] tuneR_1.4.7 brmsish_1.0.0 lunar_0.2-1 ggplot2_3.5.1
[13] ggdist_3.3.2 knitr_1.47 kableExtra_1.4.0 HDInterval_0.2.4
[17] readxl_1.4.3 posterior_1.5.0 cowplot_1.1.3 brms_2.21.0
[21] Rcpp_1.0.12 viridis_0.6.5 viridisLite_0.4.2 remotes_2.5.0
loaded via a namespace (and not attached):
[1] tensorA_0.36.2.1 rstudioapi_0.16.0 jsonlite_1.8.8
[4] magrittr_2.0.3 TH.data_1.1-2 sketchy_1.0.3
[7] estimability_1.5.1 farver_2.1.2 rmarkdown_2.27
[10] ragg_1.3.2 vctrs_0.6.5 memoise_2.0.1
[13] RCurl_1.98-1.14 htmltools_0.5.8.1 distributional_0.4.0
[16] curl_5.2.1 signal_1.8-1 cellranger_1.1.0
[19] StanHeaders_2.32.9 KernSmooth_2.23-24 htmlwidgets_1.6.4
[22] plyr_1.8.9 sandwich_3.1-0 testthat_3.2.1.1
[25] cachem_1.1.0 emmeans_1.10.2 zoo_1.8-12
[28] igraph_2.0.3 lifecycle_1.0.4 pkgconfig_2.0.3
[31] Matrix_1.7-0 R6_2.5.1 fastmap_1.2.0
[34] digest_0.6.36 colorspace_2.1-0 textshaping_0.4.0
[37] labeling_0.4.3 fansi_1.0.6 polyclip_1.10-6
[40] abind_1.4-5 compiler_4.4.1 proxy_0.4-27
[43] withr_3.0.0 backports_1.5.0 inline_0.3.19
[46] DBI_1.2.3 QuickJSR_1.2.2 pkgbuild_1.4.4
[49] highr_0.11 ggforce_0.4.2 MASS_7.3-61
[52] rjson_0.2.21 classInt_0.4-10 loo_2.7.0
[55] tools_4.4.1 units_0.8-5 ape_5.8
[58] fftw_1.0-8 glue_1.7.0 nlme_3.1-165
[61] grid_4.4.1 sf_1.0-6 checkmate_2.3.1
[64] reshape2_1.4.4 generics_0.1.3 diffobj_0.3.5
[67] gtable_0.3.5 class_7.3-22 tidyr_1.3.1
[70] tidygraph_1.3.1 xml2_1.3.6 utf8_1.2.4
[73] ggrepel_0.9.5 pillar_1.9.0 stringr_1.5.1
[76] splines_4.4.1 dplyr_1.1.4 tweenr_2.0.3
[79] lattice_0.22-6 survival_3.7-0 tidyselect_1.2.1
[82] arrayhelpers_1.1-0 gridExtra_2.3 V8_4.4.2
[85] svglite_2.1.3 stats4_4.4.1 xfun_0.45
[88] graphlayouts_1.1.1 bridgesampling_1.1-2 brio_1.1.5
[91] matrixStats_1.3.0 rstan_2.32.6 stringi_1.8.4
[94] yaml_2.3.8 boot_1.3-30 xaringanExtra_0.8.0
[97] evaluate_0.24.0 codetools_0.2-20 dtw_1.23-1
[100] ggraph_2.2.1 tibble_3.2.1 cli_3.6.3
[103] RcppParallel_5.1.7 xtable_1.8-4 systemfonts_1.1.0
[106] munsell_0.5.1 coda_0.19-4.1 svUnit_1.0.6
[109] parallel_4.4.1 rstantools_2.4.0 bayesplot_1.11.1
[112] Brobdingnag_1.2-9 bitops_1.0-7 mvtnorm_1.2-5
[115] scales_1.3.0 e1071_1.7-14 packrat_0.9.2
[118] purrr_1.0.2 crayon_1.5.3 rlang_1.1.4
[121] multcomp_1.4-25